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ABSTRACT 

We perform numerical simulations of a stellar galactic disk with initial conditions chosen to represent 
an unrelaxed population which might have been left following a merger. Stars are unevenly distributed 
in radial action angle, though the disk is axisymmetric. The velocity distribution in the simulated 
Solar neighborhood exhibits waves traveling in the direction of positive w, where u, v are the radial 
and tangential velocity components. As the system relaxes and structure wraps in phase space, the 
features seen in the ww-plane move closer together. We show that these results can be obtained also 
by a semi-analytical method. We propose that this model could provide an explanation for the high 
velocity streams seen in th e Solar neighborhood at ap proximate v in km/ s, of -60 (HR 1614), -80 
(|Arifvanto and Fuchsl[200l) . -100 (Arcturus), and -160 (jKlement et al.ll2008f) . In addition, we predict 
four new features sX v ~ -140, -120, 40 and 60 km/s. By matching the number and positions of the 
observed streams, we estimate that the Milky Way disk was strongly perturbed ~1.9 Gyr ago. This 
event could have been associated with Galactic bar formation. 
Subject headings: stellar dynamics 



1. INTRODUCTION 

The formation and evolution of galaxies is one of the 
most important topics in contemporary astrophysics. 
High-redshift cosmology provides insight into the evo- 
lution of global galaxy properties, but is unable to probe 
the internal kinematics and chemistry on sub-galactic 
scales. The Milky Way (MW), on the other hand, con- 
tains a vast amount of fossil evidence encoded in the mo- 
tions and chemical properties of its stars. The Galaxy is 
the only galaxy within which we can obtain information 
at the level of detail required to distinguish robustly be- 
tween different formation scenarios. Unfortunately, how- 
ever, the precise structure of the MW remains a topic of 
debate. In order to make progress in this field we need 
to differentiate between different Galactic models. 

As part of this effort, models aimed at explaining asym- 
metries in the Solar Neighborhood (SN) velocity space as 
the result of internal (spiral and/or bar structure) or ex- 
ternal (satellite mergers) agents have been explored in 
the past several decades. Hipparcos data r evealed a non- 
smooth local v elocity distribution of stars (|Dehnenlll998t 
IChereul et al.i ri998. 1999). These s tellar streams c annot 
simply be dissolved clusters ( Famaev et al.l 120071 ) . but 
could be caused by dynamical effects within the MW 
disk or satellite mergers. Some of these features have 
been used as tracers of non-axisymmetric Galactic disk 
structure and employed i n estimating parameters of the 
MW central ba r (Dchnen 20001: lFuxll2001: Min chev et all 
20071: IChakrabartv..2007.) and spirals dLepine el. al.ll2001t 
Quillen and Minchevl " l2005l : IChakrabartvl l2007t ). How- 
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ever, dynamical instabilities in the Galactic disk have 
been found to only relate to velocities in the range 
u,v ±50 km/s, where u,v are the radial and tangen- 
tial galactocentric velocities of SN stars, respectively. At 
this time there is no evidence that spiral or bar struc- 
ture can cause high velocity streams, thus, these are usu- 
ally attributed to merger events. For instance, the Arc- 
turus stream at v = —100 km/s, has been interpreted 
as originating from the debris of a di srupted satellite 
(jNavarro eratl l2004 fHelmi et al.l I2006D. Two recently 
disco vered streams at w ~ 80 km/s (lArifvanto and Fuch j 

I2006D and v 160 km/s ()Klement et al.l I2008D were 

assigned similar origin, based on their kinematics. There 
is plenty of evidence for past and ongoing accretion of 
small objects by the MW, the most dramatic one be- 
ing the highly dis rupted Sgr dwarf galaxy identified by 
(|Ibata et al.lll994 fl995). But what about clues of more 
massive MW mergers? 

Cosmological simulations show that massive minor 
mergers are likely to have happened during the lifetime of 
a MW-sized host system. Early generations of nurn erical 
simulations (jOuinn "eTall [19931 : 1 Walker et al.l[l996h have 
investigated the heating of galactic disks by mergers, 
finding that disks heated in this way are similar to the 
MW thick disk. By considering realistic satellite orbits 
and conditions, more recent attempts to model disk heat- 
ing in a cosmological context suggest sm aller efficiency 
than previously thought. For instance, iHopkins et al.l 
1^08) estimated that the Milky Way could have survived 
as many as ^5-10 minor (mass ratio to host disk ^1:10) 
mergers in the last 10 Gyr. By estimating the maximum 
mass ratio merger for a MW-sized galaxy, the authors 
point out that the Galaxy could even have survived merg- 
ers o f mass ratio ^ 1 : 4 — 1 : 3 without destroying the 
disk. iKazantzidis et aO ()2008l ) estimated the number of 
massive subhaloes accreted since z~ 1 to be at least one 
object with a mass ~ Mdisc and five objects more massive 
than 20% Mdisc- These findings are consistent with simu- 
lations by other groups {B enson et al. 2004; Stochr 200d 
iDe Lucia and Helmil [20081 : IVillalobos and Helmil l2008l : 



2 



IStewart et all 120081 : iPurcell et all I2009D . Even though 
massive minor mergers may be less frequent, they are 
able to reach the center of the host system (provided 
they are also dense) thanks to dynamical friction, caus- 
ing important changes in the structure and kinematics 
of the host disc. By usi ng data from the Two Mic ron 
All-Sky Survey r2MASS'l. ICole and WeinbeT^ (|200l es- 
timated that the Milky Way bar is likely to have formed 
more recently than 3 Gyr ago and suggested that this 
event could have been triggered by a now-merged satel- 
lite. Formation of central bars as the result of satellite- 
disk en counters has also been observed in N-body s imu- 
lations (jWalker et al.lll996HKa"zantzidis et al.ll2008l ). 

All this evidence from both observations and simu- 
lations implies that merging satellites could contribute 
to the heating of the MW disc. Is it possible to re- 
late the effects of such an event to features in the 
local velocity space (uv-plane)? Instead of interpret- 
ing streams as debris from disrupted small sa t ellites , 
as done by [Na varro ct al.' (2004); Hclmi ct al. (" 20061 ): 
I Arifvanto andFuc hs (2006); Klcmcnt ct al. (200^, we 
ask a different question here: Is it possible that some 
overdensities in the SN velocity space are simply the re- 
sponse of the MW disk to the sudden energy kick im- 
posed by a massive satellite in the past? 

The signature of this perturbation will be present in 
the MW stellar kinematics during the relaxation time 
after the event. Depending on the time of this event, it 
is possible that the Galactic disk is still undergoing relax- 
ation, although at this time the impact imprint may be 
extremely faint. However, because stars of the thick disk 
spend relatively little time near the galactic plane, where 
the spiral arm heating and scattering by giant molecular 
clouds is most vigorous, radial mixing within the thick 
disk is unlikely to wipe out the signature of a past event 
too quickly. 

In this letter we look at the evolution of a non-relaxed 
galactic disk and its manifestation in the velocity distri- 
bution of a simulated SN. 

2. SIMULATIONS OF AN UNRELAXED DISK 

At present it is not computationally feasible to achieve 
the statistics needed for resolving a small spatial region 
with N-body simulations, even in the case of a single 
galaxy. Thus, a realistic galaxy collision simulation does 
not provide the resolution needed for our purposes, i.e., 
resolve a fictitious SN to look for structure in the local 
velocity space. We choose initial velocities for particles 
axisymmetrically by means of Gaussian distributions in 
u and V with corresponding standard deviations cr„ and 
(7^, consistent with a hot stellar population. However, we 
purposely choose them so that they are not evenly dis- 
tributed in their radial oscillation. This serves as a proxy 
for choosing a population that is unrelaxed or unevenly 
distributed in phase space, such as might be left after 
a merger. We further simplify the problem by consider- 
ing only two dimensions, assuming the vertical motion 
of stars is decoupled from the motion in the plane of 
the Galaxy. We are mainly concerned with a flat rota- 
tion curve but also discuss the effect of a decreasing and 
increasing one in section 21 In addition to an axisymmet- 
ric system we also simulate a disk perturbed by a cen- 
tral bar as a pure quadrupole. A detailed description o f 
the perturbation can be found in lMinchev et al.l (|2007l ) . 



To explore the time development of the system, we do 
not time-average over position and velocity ve ctors, as 
it is frequently done in test-partic le simulations (|Dehnenl 
2000; Fux 2001- |Minchev and Qui llcn 2007) where no dy- 
namical development of the system is expected. To con- 
vert to real units we use Local Standard of Rest (LSR) 
tangential velocity of 220 km/s, and Galactocentric dis- 
tance of 7.8 kpc. 

Note that the Gaussian distributions in u and v pro- 
vide initial conditions (ICs) sampled non-uniformly in 
9, the radial epicyclic angle. These ICs were found to 
induce an initial radial expansi on in the disc consistent 
with the N-b o dy sim ulations bv lQuinn et al.l (|1993l ) and 
I Walker et al.l (|1996D . where they found that the host 
disks respond by spreading both radially and vertically. 
However, after a couple of rotations the density distribu- 
tion appears smooth and is axisymmetric. 

3. RESULTS 
3.1. Density waves in velocity space 

We now look for the effect of our ICs on a simulated 
SN velocity distribution. In figure [1] we show the time 
development of the local velocity field in three differ- 
ent ways. The first row shows the ww-plane (contours) 
and the v distribution (solid line). The second row plots 
it^ -|- 2v'^y /^ versus v as done by lArifvanto and Fuchsl 
2006| ) and iKlement et~al\ ()2008( ). In this simulation 
cr„ — 50 km/s. The sample shown is limited to a ra- 
dius of 100 pes around our fictitious Sun. We show six 
time outputs up to t = 10 rotations at tq. Note that 
features in the ww-plane get closer together as time in- 
creases. We interpret this as wrapping in phase space 
on a timescale associated with the epicyclic frequency. 
The features are not ori ented along constan t eccentricity 
surfaces as predicted bv lHelmi et al.l (|2006D for particles 
trapped from a minor merging galaxy. In our case the 
arcs are oriented in the opposite direction, since they are 
centered on {u, v — vq) = 0, in other words these are 
constant energy surfaces. 

Note that the features in the ww-plane are manifested 
in the tail of the w-distribution, as well as in (u^-|-2'i;^)^/^, 
as overdensities traveling in the direction of positive v. 
In section 13.31 these are used to match to high velocity 
streams observed in the SN stellar population. 

It is important to make a point here concerning the 
average distance of our sample from the Sun. Because 
of the differential rotation of the Galaxy, at radii inte- 
rior to ro the relaxation is completed faster, whereas the 
opposite is true for stellar samples at r > tq. Thus, for 
a given time, the separation of features in the uw-plane 
depends on the Galactocentric distance of our sample. 
Consequently, as sample depth increases we sample a 
large range of Galactic radii causing the waves to inter- 
fere and either enhance or (mostly) wipe out structure 
in the ww-plane. 

To better understand the mechanism giving rise to 
these ripples in velocity space, next we use the epicyclic 
approximation to reproduce the numerical results just 
described. 

3.2. Semi- analytical approach 
We assume a flat rotation curve with the parameter 
7 = 2D, /k — \/2, where f2 and k are the angular rotation 
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Fig. 1. — Time development of the local velocity field for an axisymmetric disk with ICs as described in the text, presented in 
three different ways. First row shows the uv-plane (contours) and the tangential velocity distribution (solid line). The second row plots 
(u^ + 2v'^)^^'^ versus v. The sample shown is limited to a radius of 100 pc from our fictitious Sun. We show six time outputs up to t = 10 
SN rotations. As time increases features get closer together as phase wrapping takes place. 
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Fig. 2. — The uv-plane computed using a weighting function utilizing eq. 30b from IDehnenI l|1999l ). Note the striking similarity to our 
numerical simulations (figure [TJ. 



rate and epicyclic frequency of a particle in a circular or- 
bit. Since we work in units of and tq, the gravitational 
potential for a flat rotation curve <I>(7-) — v'^ln{r/ro) is 
zero for stars in the SN. For every position on the uv- 
plane energy and angular momentum can be computed 
a s E{u, v) = ^{u ^ + (1 + v)"^) and L(v) = 1 + i>. 

IDehnenI (|1999|) derived a good approximation for the 
frequency of epicycles accurate at large epicyclic ampli- 
tude (his equation 30d): 

where 

= 1- L^/L^{E) (2) 

is the square of the orbital eccentricity (his equation 30b). 
These functions depend on k{E) and L{E), which are 
the epicyclic frequency and the angular momentum for a 
circular orbit with energy E 



«(£;) = %/2cxp ii^-E 



L(£')==cxp 



(3) 
(4) 



where we have used expressions from Table 1 bv lDehiienI 
([1999), given in units of the circular velocity. Since lor 
depends on L and E, it can be computed for every posi- 
tion on the wt;-plane in the SN. 



We assume that our initial particle distribution is not 
evenly distributed in the angle associated with epicyclic 
motion, 9. If the initial phase space density distribution 
is skewed along this angle at 6*0, then at a time At later 
there will be a maximum of particles at 



e{L,E) 



-UB.iL,E)At 



(5) 



To mimic the effect of this we construct a weighting func- 
tion that gives a maximum at 0{L, E) = mod 27r. 



w{L,E) = ^{l 



COs{ujRt)) 



(6) 



Our initial phase space particle distribution is more un- 
evenly distributed in action angle at large eccentricity so 
we can improve our weighting function by multiplying 
the cosine function by the square of the orbital eccen- 
tricity. 



;(L, E) = exp(-e/eo)(l + cos{ujRt)), 



(7) 



where the exponential function mimics a Gaussian veloc- 
ity dispersion and cq is a constant. In the above equation 
lor and e depend on L,E. Equation [7] can be computed 
for every position on the uw-plane for different values of 
t. The result is shown in figure [2] and looks remarkably 
close to what is seen in the test particle simulations. 

3.3. Relating to moving groups 
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We now discuss the possibility that ringing in the 
Galaxy is the reason for four high velocity streams ob- 
served in the local velocity field. We search for a par- 
ticular time in our synthetic velocity distributions, for 
which we get a satisfactory match to all four. Figure 
[3] presents velocity field plots from our numerical sim- 
ulations for time t = 8.67 rotations at tq, in three dif- 
ferent ways as in figure [T] The dashed lines indicate 
the position of the observed overdensities. From left to 
right these are as follows: (1) A high velocity stream 
at V = —16 lb 20 km/s (h e reafte r STRl) was recently 
reported by iKlement et al.l (|2008D which, based on its 
kinematics, is thought to belong to the thick disk. (2) 
Arcturus is a moving group lagging the LSR by 100 km/s. 
Its metal-poor nature and significant age are consistent 
wi th the thick disk. A d etailed investigation of its origin 
bv I Williams et all (|2008f ) found that the chemical results 
are consistent with a dynamical origin but do not entirely 
rule out a merger one. T he upper left-hand panel in fig- 
ure 1 bv lWilliams et all (j2008l) presents a uw-plot of the 
Arcturus stream. Centered on a narrow v, it spreads over 
the range —100 < u < 100, consistent with our results 
(left panel in figure [3]). (3) A stream with characteristics 
appropriate for the thick di sk atv = —80 km/s (he reafter 
STR2) was found by lArifvanto and Fuchsl (|2006[ ) using 
data extracted from various catalogues. (4) The moving 
group HR 1614 at v = -60 km/s, is thought to be a dis- 
persed ogen cluster because of its chemical homogeneity 
(|Eggenlll99^ iDe Silva et aIll2007D at a distance of 40 pc 
from the sun. It is intriguing that, similarly to Arcturus, 
in the uv-planc this stream has a spread in u, again con- 
sistent wit h figure Bl However , —50 < u < 20 km/s 
(figure 5 in IDe Silva et al.|[2007f ). i.e., the corresponding 
wave giving rise to HR 1614 is distorted toward negative 
u. This is consistent with the effect of the bar, given 
the proximity of this overdensity to the Hercules stream. 
Note that our interpretation for HR 1614's location at 
V = 60 km/s does not contradict the possibility of it be- 
ing a dispersed cluster, as long as it is older than the 
time of the merger event, since the streams in the model 
proposed here are only defined by their kinematics. 

In addition to the four observed streams, our new 
model predicts overdensities at approximately every 20 
km/s. However, in the range —50 < w < 50 km/s, 
the central peak of the velocity distribution dominates. 
Thus, we are left with four new, easily identifiable over- 
densities at V —140, -120, 40 and 60 km/s, indicated 
by the dotted lines in figure O 

In order to look for the streams our r nodel predicts, we 
combi ned t he observational samp les bv lNordstrom et al.l 
(|2004[ ) and lSchuster et all (|2006ft . We select a thick disk 
dominated sample by considering the metallicity range 
— 1.1 <[Fe/H]< —0.55 dex. In figure [4] we present the 
result for sample depths d,nax = 80 and dmax = 150 
pc in the same fashion as the left-hand panel of fig- 
ure [S] With this small number of stars (N=451,766 for 
dmax — 80, 150 pc), this result is not statistically signif- 
icant to provide convincing evidence for the validity of 
our model. However, the resemblance of the observed 
w-distribution to the ones resulting from our simulations 
and semi-analytical approach, is striking. 



4. DISCUSSION AND CONCLUSIONS 
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Fig. 3. — Three ways of plotting the same data, as done in 
figure [T] for a particular time matching four observed high velocity 
local streams (dashed lines). From left to right, in km/s, these 
are STRl at d Ri -160, Arcturus a.t v -100, STR2 at v ^ -80 
and HR 1614 at v fa -60. The dotted lines at u -140, -120, 
40 and 60 km/s show the positions of four new stream this new 
model predicts. Note that the density oscillations giving rise to 
these features have an amplitude of less than one percent. 
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Fig. 4. — O b servat ional data form lNordstrom et al. > f2004') and 
ISchuster et al.l 120061 ) combined samples. We select stars in the 
metallicity range —1.1 <[Fe/H]< —0.55 dex. The dashed and dot- 
ted lines indicate the observed and predicted overdensities as in 
figure \3\ 



We have shown that an axisymmetric galactic disk sub- 
jected to an initial energy kick approximating a massive 
Galactic merger induces waves in the SN velocity field 
propagating in the direction of positive v, which appear 
as overdensities in the tail of the tangential velocity dis- 
tribution. By comparing our synthetic models to obser- 
vations, a satisfactory match to four SN high velocity 
streams is achieved toward the end of the disk relax- 
ation at t — 8.67 SN rotations (figure [3]). In addition, 
we predict the existence of four (or more) new features 
at V fa -140, -120, 40 and 60 km/s. Our results allow 
us to make an estimate for the time of the event. For a 
Galactocentric distance of 7.8 kpc and a LSR tangential 
velocity of 220 km/s, t — 8.67 rotations corresponds to 
~ 1.9 Gyr. If our model is correct, then all observed and 
predicted streams must be older than the time of the disk 
stirring event. All four of the known features discussed in 
section [331 have ages >~ 2 Gyr, which is consistent with 
our prediction for the time of the impact. This model 
also arg ues against purely diffusi ve stochastic heating 
models (iJenkins and Binncv 1990l: ISellwood and BinnevI 
I2002t iMinchev and'Quillen,2006i) for the Galactic disk. 

In addition to a flat rotation curve, we have also con- 
sidered a power law initial tangential velocity = 
vo{r/ro)^ with (3 — 0.2,-0.2 corresponding to a rising 
and a declining rotation curve, respectively. We found 
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that the separation of the features in the uu-plane de- 
creases with time more rapidly for /3 — —0.2 and more 
slowly in the case of /3 = 0.2, as expected. However, at rp 
our results remained the same. The ^ 20 km/s separa- 
tion of features in the local velocity field arises naturally 
as a result of the Galactocentric distance and tangential 
velocity of the LSR, i.e., it is determined only by the 
LSR angular velocity. Thus, not only is this model inde- 
pendent of the MW rotation curve, but it can be used to 
provide constraints on /? by observations of the velocity 
field at Galactic radii different than rg. 

Features in the uv-plane represent curves of constant 
energy and are oriented opposite t o the constant eccen- 
tricity curves in iHelmi et al.1 (|2006D 's model. As more SN 
stars are surveyed with RAVE and GAIA, the shape of 
these will be better resolved and we will be able to tell the 
difference between Helmi's model and ours. We predict 
a shift in location of features as a function of Galacto- 
centric radius (closer together for shorter radii) as the 
distance between the features depends on the epicyclic 
frequency. Deeper surveys, such as ARGOS, SEGUE, 
BRAVA, and AP OGEE, could search fo r this s hift in the 
way sugg ested bv lMinchev and QuillenI (|2008l ). 

We have checked that the growth, or longer term ef- 
fect, of a central bar does not cause similar features in 
the wii-plane, thus the bar is not responsible for such ra- 
dial perturbations. However, an increase in central mass 
associated with a merger could cause such variations in 
the epicyclic action-angle distribution. From recent cos- 
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mological simulations it is now known that satellites of 
mass ratio to host disk > 1 : 10 merge on highly ec- 
centric , nearly radial orbit s in a couple of dynamical 
times (jHopkins et al.|[2008f) . Then they quickly merge 
by dumping their mass in the center of the host disk. It 
would appear to be simple to look for initial conditions 
from an N-body simulation. However this is actually 
nontrivial for a number of reasons: (1) not enough parti- 
cles to resolve high velocity structure; (2) large range of 
possibilities; (3) time dependent phenomena associated 
with mergers. 

Our short timescale is consistent with the esti- 
mated age of th e Gal actic bar, as measured by 
iCole and Weinberg! (l200l) (< 3 Gyr ago). This suggests 
that the same event that caused the formation of the 
Galactic bar could have left the stellar disk unrelaxed, 
thus giving rise to the observed high velocity streams. 
Another po ssibility for stirring up the MW disk is th e 
wCen event (jBekki and Freemanll2003HMeza et al.l2005H . 

Our choice of ICs for an unrelaxed disk is arbitrary. 
N-body simulations could be explored to see if such a 
distribution could arise from a merger and motivate bet- 
ter choices of ICs for future particle integration studies. 
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